##Title:Main Analysis
##Date: 18 Sept 2023

```{r preparation}
library(tidyverse)

library(plm)
library(modelsummary)

dat <- read.csv("D:/Essex_PhD/3rd_year/Data for R&R/Replication_R1_20240918/replication_20231015.csv") 


```



```{r Table 1 Distribution of Observations of Grids with Peacekeeping and/or Aid}
dat <- dat%>%
  mutate(catergory = case_when(
    treat_aid ==1 & treat_pko == 1 ~ "3",
    treat_aid ==0 & treat_pko == 1 ~ "2",
    treat_aid ==1 & treat_pko == 0 ~ "1",
    treat_aid ==0 & treat_pko == 0 ~ "0"
  ))

table(dat$catergory)

table(dat$osv)
#All observation -#OSv with 0 = #GRID with OSV
#31908-30880 = 1028
```


```{r  table2a log_OSV(t+1) rescale some variables}
#All sample
dat$prec_gpcp_100 <- (dat$prec_gpcp)/100
dat$conflict_intensity_100 <- (dat$conflict_intensity)/100
dv_c_iv_b.1 <- plm(log_osv_forward1 ~ treat_pko + treat_aid + 
                        treat_pko*treat_aid, 
                         data = dat,
                         index = c("gid", "year"), model = "within")

dv_c_iv_b.2 <- plm(log_osv_forward1 ~ treat_pko + treat_aid + 
                     treat_pko*treat_aid +
                     lnpop + prec_gpcp_100 + conflict_intensity_100 + 
                     spatial_lag_osv + spatial_lag_pko +
                     task_tamm + humaid_assist_tamm,
                     data = dat,
                         index = c("gid", "year"), model = "within")

# Reb sample

reb_dv_c_iv_b.1 <- plm( log_reb_osv_forward1 ~ treat_pko + treat_aid + 
                 treat_pko*treat_aid , 
                         data = dat,
                         index = c("gid", "year"), model = "within")

reb_dv_c_iv_b.2 <- plm( log_reb_osv_forward1 ~ treat_pko + treat_aid + 
                 treat_pko*treat_aid + 
                 lnpop  + prec_gpcp_100 + conflict_intensity_100 + 
                 spatial_lag_reb_osv + spatial_lag_pko + 
                     task_tamm + humaid_assist_tamm, 
                         data = dat,
                         index = c("gid", "year"), model = "within")

# GOV sample

gov_dv_c_iv_b.1 <- plm(log_gov_osv_forward1 ~ treat_pko + treat_aid + 
                 treat_pko*treat_aid , 
                         data = dat,
                         index = c("gid", "year"), model = "within")

gov_dv_c_iv_b.2 <- plm( log_gov_osv_forward1 ~ treat_pko + treat_aid + 
                 treat_pko*treat_aid + 
                 lnpop  + prec_gpcp_100 + conflict_intensity_100 + 
                 spatial_lag_gov_osv + spatial_lag_pko+ 
                     task_tamm + humaid_assist_tamm , 
                         data = dat,
                         index = c("gid", "year"), model = "within")



modelsummary(list( dv_c_iv_b.1, dv_c_iv_b.2,
                   reb_dv_c_iv_b.1, reb_dv_c_iv_b.2,
                   gov_dv_c_iv_b.1, gov_dv_c_iv_b.2), 
             stars = TRUE)
```




```{r table2b OSV(t+1) Binary rescale some variable}
#All sample
dv_b_iv_b.1 <- plm(osv_binary_forward1 ~ treat_pko + treat_aid + 
                        treat_pko*treat_aid, 
                         data = dat,
                         index = c("gid", "year"), model = "within")

dv_b_iv_b.2 <- plm(osv_binary_forward1 ~ treat_pko + treat_aid + 
                        treat_pko*treat_aid +
                     lnpop + prec_gpcp_100 + conflict_intensity_100 + 
                     spatial_lag_osv + spatial_lag_pko + 
                     task_tamm + humaid_assist_tamm,
                     data = dat,
                         index = c("gid", "year"), model = "within")

#Reb sample

reb_dv_b_iv_b.1 <- plm(reb_osv_binary_forward1 ~ treat_pko + treat_aid + 
                 treat_pko*treat_aid , 
                         data = dat,
                         index = c("gid", "year"), model = "within")

reb_dv_b_iv_b.2 <- plm(reb_osv_binary_forward1 ~ treat_pko + treat_aid + 
                 treat_pko*treat_aid + 
                 lnpop + prec_gpcp_100 + conflict_intensity_100 + 
                 spatial_lag_reb_osv + spatial_lag_pko + 
                     task_tamm + humaid_assist_tamm, 
                         data = dat,
                         index = c("gid", "year"), model = "within")

#Gov sample
gov_dv_b_iv_b.1 <- plm(gov_osv_binary_forward1 ~ treat_pko + treat_aid + 
                 treat_pko*treat_aid , 
                         data = dat,
                         index = c("gid", "year"), model = "within")

gov_dv_b_iv_b.2 <- plm( gov_osv_binary_forward1 ~ treat_pko + treat_aid + 
                 treat_pko*treat_aid + 
                 lnpop  + prec_gpcp_100 + conflict_intensity_100 + 
                 spatial_lag_gov_osv + spatial_lag_pko+ 
                     task_tamm + humaid_assist_tamm , 
                         data = dat,
                         index = c("gid", "year"), model = "within")

modelsummary(list( dv_b_iv_b.1, dv_b_iv_b.2,
                   reb_dv_b_iv_b.1, reb_dv_b_iv_b.2,
                   gov_dv_b_iv_b.1, gov_dv_b_iv_b.2), 
             stars = TRUE)
```



```{r Figure 2}
gid.fact<-as.factor(dat$gid)
year.fact<-as.factor(dat$year)
treat_aid.fact <- as.factor(dat$treat_aid)
PKO_presence <- as.factor(dat$treat_pko)

model_same <- lm(osv_binary_forward1 ~ treat_aid.fact + PKO_presence + treat_aid.fact:PKO_presence + 
                 lnpop + prec_gpcp + conflict_intensity + 
                 spatial_lag_osv + spatial_lag_pko + 
                 gid.fact + year.fact + 
                     task_tamm + humaid_assist_tamm,
                 data = dat)

library(effects)

inter.plot1 <- plot(Effect(c("PKO_presence", "treat_aid.fact"),
                          mod= model_same), 
     lines= list(multiline=TRUE,  lty= 2:3),
     confint=list(style="auto"),
     xlab = "PKO presence (binary)",
     ylab =  "One-sided violence",
     main = "Interaction plot of Aid*PKO presence (binary) plot")
```

```{r Figure 3}
#Figure rebel

reb_model_same <- lm(reb_osv_binary_forward1 ~ treat_aid.fact + PKO_presence + treat_aid.fact:PKO_presence + 
                 lnpop + prec_gpcp + conflict_intensity + 
                 spatial_lag_osv + spatial_lag_pko + 
                 gid.fact + year.fact + 
                     task_tamm + humaid_assist_tamm,
                 data = dat)

library(effects)
reb_inter.plot1 <- plot(Effect(c("PKO_presence", "treat_aid.fact"),
                          mod= reb_model_same), 
     lines= list(multiline=TRUE,  lty= 2:3), confint=list(style="auto"),
     xlab = "PKO presence (binary)",
     ylab =  "Reb One-sided violence")

#Figure Gov
gov_model_same <- lm(gov_osv_binary_forward1 ~ treat_aid.fact + PKO_presence + treat_aid.fact:PKO_presence + 
                       lnpop + prec_gpcp + conflict_intensity + 
                       spatial_lag_osv + spatial_lag_pko + 
                       gid.fact + year.fact + 
                       task_tamm + humaid_assist_tamm,
                     data = dat)

gov_inter.plot1 <- plot(Effect(c("PKO_presence", "treat_aid.fact"),
                               mod= gov_model_same), 
                        lines= list(multiline=TRUE,  lty= 2:3), confint=list(style="auto"),
                        xlab = "PKO presence (binary)",
                        ylab =  "Gov One-sided violence")

```

```{r Table 3}

#DV log

dv_c_iv_non <- plm(log_osv_forward1 ~ catergory + 
                     lnpop + prec_gpcp + conflict_intensity + 
                     spatial_lag_osv + spatial_lag_pko + 
                     task_tamm + humaid_assist_tamm,
                     data = dat,
                         index = c("gid", "year"), model = "within")

reb_dv_c_iv_non <- plm(log_reb_osv_forward1~ catergory + 
                     lnpop + prec_gpcp + conflict_intensity + 
                     spatial_lag_osv + spatial_lag_pko + 
                     task_tamm + humaid_assist_tamm,
                     data = dat,
                         index = c("gid", "year"), model = "within")

gov_dv_c_iv_non <- plm(log_gov_osv_forward1 ~ catergory + 
                     lnpop + prec_gpcp + conflict_intensity + 
                     spatial_lag_osv + spatial_lag_pko + 
                     task_tamm + humaid_assist_tamm,
                     data = dat,
                         index = c("gid", "year"), model = "within")


#DV binary

dv_b_iv_non <- plm(osv_binary_forward1 ~ catergory + 
                     lnpop + prec_gpcp + conflict_intensity + 
                     spatial_lag_osv + spatial_lag_pko + 
                     task_tamm + humaid_assist_tamm,
                     data = dat,
                         index = c("gid", "year"), model = "within")

reb_dv_b_iv_non <- plm(reb_osv_binary_forward1 ~ catergory + 
                     lnpop + prec_gpcp + conflict_intensity + 
                     spatial_lag_osv + spatial_lag_pko + 
                     task_tamm + humaid_assist_tamm,
                     data = dat,
                         index = c("gid", "year"), model = "within")

gov_dv_b_iv_non <- plm(gov_osv_binary_forward1 ~ catergory + 
                     lnpop + prec_gpcp + conflict_intensity + 
                     spatial_lag_osv + spatial_lag_pko + 
                     task_tamm + humaid_assist_tamm,
                     data = dat,
                         index = c("gid", "year"), model = "within")

modelsummary(list(dv_c_iv_non, reb_dv_c_iv_non, gov_dv_c_iv_non,
                  dv_b_iv_non, reb_dv_b_iv_non, gov_dv_b_iv_non),
             stars = TRUE)
```

```{r Figure 4}
library(sensemakr)
dat$pko_aid <-(dat$treat_aid)*(dat$treat_pko)


#all sample

#FE with pko_aid
sens_fe.1 <- lm(osv_binary_forward1 ~ treat_aid + treat_pko + pko_aid + 
                                   lnpop + prec_gpcp  + conflict_intensity + 
                                   spatial_lag_osv + spatial_lag_pko + 
                                    task_tamm + humaid_assist_tamm +
                                   as.factor(gid) + as.factor(year),   data = dat)
#FE without pko_aid
sens_fe.2 <- lm(osv_binary_forward1 ~ treat_aid + treat_pko + 
                                   lnpop + prec_gpcp  + conflict_intensity + 
                                   spatial_lag_osv + spatial_lag_pko + 
                                   task_tamm + humaid_assist_tamm +
                                   as.factor(gid) + as.factor(year),   data = dat)

sens.fe.1 <- sensemakr(model = sens_fe.1, 
                         treatment = "pko_aid",
                         benchmark_covariates = c( "conflict_intensity"),
                         kd = 1:3)

sens.fe.2 <- sensemakr(model = sens_fe.2, 
                         treatment = "treat_pko",
                         benchmark_covariates = c( "conflict_intensity", "treat_aid"),
                         kd = 1:3)
par(mfrow=c(1,2))
plot(sens.fe.1, lim = 0.05, lim.y = 0.05)
plot(sens.fe.2, lim = 0.05, lim.y = 0.05)

#reb sample

#FE with pko_aid
reb_sens_fe.1 <- lm(reb_osv_binary_forward1 ~ treat_aid + treat_pko + pko_aid + 
                                   lnpop + prec_gpcp  + conflict_intensity + 
                                   spatial_lag_osv + spatial_lag_pko + 
                                    task_tamm + humaid_assist_tamm +
                                   as.factor(gid) + as.factor(year),   data = dat)
#FE without pko_aid
reb_sens_fe.2 <- lm(reb_osv_binary_forward1 ~ treat_aid + treat_pko + 
                                   lnpop + prec_gpcp  + conflict_intensity + 
                                   spatial_lag_osv + spatial_lag_pko + 
                                   task_tamm + humaid_assist_tamm +
                                   as.factor(gid) + as.factor(year),   data = dat)

reb_sens.fe.1 <- sensemakr(model = reb_sens_fe.1, 
                         treatment = "pko_aid",
                         benchmark_covariates = c( "conflict_intensity"),
                         kd = 1:3)

reb_sens.fe.2 <- sensemakr(model = reb_sens_fe.2, 
                         treatment = "treat_pko",
                         benchmark_covariates = c( "conflict_intensity", "treat_aid"),
                         kd = 1:3)
par(mfrow=c(1,2))
plot(reb_sens.fe.1, lim = 0.05, lim.y = 0.05)
plot(reb_sens.fe.2, lim = 0.05, lim.y = 0.05)


#gov sample

#FE with pko_aid
gov_sens_fe.1 <- lm(gov_osv_binary_forward1 ~ treat_aid + treat_pko + pko_aid + 
                      lnpop + prec_gpcp  + conflict_intensity + 
                      spatial_lag_osv + spatial_lag_pko + 
                      task_tamm + humaid_assist_tamm +
                      as.factor(gid) + as.factor(year),   data = dat)
#FE without pko_aid
gov_sens_fe.2 <- lm(gov_osv_binary_forward1 ~ treat_aid + treat_pko + 
                      lnpop + prec_gpcp  + conflict_intensity + 
                      spatial_lag_osv + spatial_lag_pko + 
                      task_tamm + humaid_assist_tamm +
                      as.factor(gid) + as.factor(year),   data = dat)

gov_sens.fe.1 <- sensemakr(model = gov_sens_fe.1, 
                           treatment = "pko_aid",
                           benchmark_covariates = c( "conflict_intensity"),
                           kd = 1:3)

gov_sens.fe.2 <- sensemakr(model = gov_sens_fe.2, 
                           treatment = "treat_pko",
                           benchmark_covariates = c( "conflict_intensity", "treat_aid"),
                           kd = 1:3)
par(mfrow=c(1,2))
plot(gov_sens.fe.1, lim = 0.05, lim.y = 0.05)
plot(gov_sens.fe.2, lim = 0.05, lim.y = 0.05)
```

